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Abstract 

We present numerical results for the location of the chiral critical line at finite 
temperature and zero and non-zero baryon density for QCD with Nf = 2 + 1 
flavours of staggered fermions on lattices with temporal extent N t = 4. For 
degenerate quark masses, we compare our results obtained with the exact RHMC 
algorithm with earlier, inexact R-algorithm results and find a reduction of 25% 

Q-i, 

in the critical quark mass, for which the first order phase transition changes to 
a smooth crossover. Extending our analysis to non-degenerate quark masses, we 
map out the chiral critical line up to the neighbourhood of the physical point, 
^ ■ which we confirm to be in the crossover region. Our data are consistent with a 

tricritical point at (m u ^ = 0, m s ~500) MeV. 

We also investigate the shift of the critical line with finite baryon density, 
by simulating with an imaginary chemical potential for which there is no sign 
problem. We observe this shift to be very small or, conversely, the critical end- 
point n c (m u ^d,m s ) to be extremely quark mass sensitive. Moreover, the sign of 
this shift is opposite to standard expectations. If confirmed on a finer lattice, 
it implies the absence of a critical endpoint for physical QCD at small chemical 
potential, or another revision of the QCD phase diagram. We critically examine 
earlier lattice determinations of the QCD critical point, and find them to be 
in no contradiction with our conclusion. Hence we argue that finer lattices are 
required to settle even the qualitative features of the QCD phase diagram. 



1 Introduction 



Based on the property of asymptotic freedom, a fundamental prediction of QCD with 
three flavours of quarks is the transition from the familiar hadronic physics at low 
temperatures to a regime of "deconfmed" quark gluon plasma at high temperatures. 
Whether this transition is characterised by singular behaviour of the partition func- 
tion corresponding to a first or second order phase transition, or merely represents a 
smooth and analytic crossover between different dynamical regimes, depends crucially 
on the choice of the quark masses and the net baryon density specified by its chemical 
potential, hb- In the following we shall assume the light quarks to be degenerate, 
m u = rrid = m B) d, and vary m s independently. We couple the quark chemical potential 
\i = He 1 3 to the light quarks only, except for the degenerate case Nf = 3, where 
all quarks are coupled. The parameter space of the theory considered here is thus 
four- dimensional, {m u ^, m s , T, /ib}. 

The first task in determining the phase diagram in this parameter space consists of 
finding the (pseudo-)critical temperature T (m u ^ d ,m s , fi B ), defined e.g. by the peak of 
some susceptibility, which represents the boundary between the hadronic and plasma 
regimes. The independent variables m u ^, m s , fiB then span a 3d parameter space with 
regions of first order phase transitions and analytic crossover separated by surfaces 
of second order phase transitions. In order to identify the order of phase transitions, 
and the location of the critical surfaces in particular, finite size scaling analyses are 
necessary. 

Let us first discuss the situation for = 0, shown schematically in Fig. [T] (for early 
references, see, e.g. [T]). Gauge invariant, local order parameters characterising the 
transition only exist in the extreme cases of zero or infinite quark masses, namely the 
chiral condensate and the Polyakov loop, respectively. These limiting theories thus 
must feature singular phase transitions, and one may write down effective theories 
of the Ginzburg-Landau type for the order parameters [2J. It is numerically well- 
established that the phase transition is first order in the quenched jl] limit, and there 
is strong numerical evidence for first order in the chiral [3] limit. Since first order 
phase transitions are robust against small variations of the parameters of the theory, 
the first order regions must extend by a finite amount into the quark mass plane. On 
the other hand, simulations have revealed smooth crossover behaviour for intermediate 
quark masses, which implies second order boundary lines between the first order and 
crossover regions. 

In the case of heavy dynamical quarks, the relevant symmetry is the Z3 center 
symmetry, and the weakening, by the dynamical quarks, of the first-order transition 
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Figure 1: Schematic phase transition behaviour of three flavour QCD for different 
choices of quark masses (from 111), at zero density. 

which occurs in the Yang-Mills theory, is understood qualitatively [5] and to a large 
extent quantitatively. Simulations have determined the second-order line to correspond 
to a meson mass of about 2 GeV [5J, and have confirmed the expectation that the 
universality class is that of the 3d Ising model. In the case of light quarks, numerical 
simulations are more difficult, and very little is known quantitatively about the location 
of the second order boundary line. The only point computed to some accuracy with 
standard staggered fermions is the chiral critical point 0m M ^ = m s = m c (fi = 0) = 
on the Nf = 3 diagonal [H El [H] , which also was determined to belong to the 3d Ising 
universality class [7]. 

While the statement about the universality class concerns infrared physics and thus is 
stable against cut-off effects, the location of the critical point in the physical mass plane 
turns out to be very strongly affected. To date calculations have been performed on 
lattices with 4 time-slices only (N t = 4, implying a lattice spacing a = ~ 0.3 fm), 
but simulations with improved actions indicate values for 777.Q, and the associated pion 
mass, which are considerably smaller than the standard action result [7j. Moreover, 
all these simulations used the so-called R-algorithm [10] , which has stepsize errors and 
therefore gives only approximate results in the absence of a careful extrapolation to 
zero stepsize. In any case, all current results are consistent with the physical point 
being in the crossover regime. 

In the presence of a chemical potential the second order boundary lines turn into sur- 
faces, as indicated in Fig. [2j The qualitative features of the (T — /i) phase diagram now 
depend crucially on the curvature at fi = 0, d 2 m c /dfi 2 (0). The common expectation is 
that this curvature is positive. Hence the physical point, once the chemical potential is 
increased, will be closer to the critical line, and intersect it for a critical chemical po- 
tential /i c . For values larger than /z c a first order phase transition is expected. Clearly, 
this is not the case for negative curvature of the critical surface. 

1 The superscript "c" here and in the following refers to "critical", not to charm. 
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Figure 2: The chiral critical surface in the case of positive and negative curvature. If 
the physical point is in the crossover region for \i = 0, a finite \i phase transition will 
only arise in the scenario with positive curvature. 

In this work we present a comprehensive numerical study mapping out the chiral 
critical line in simulations of the standard staggered action on several lattices with 
N t = 4. Upon repeating the computation for the Nf = 3 chiral critical point with 
the rational hybrid Monte Carlo (RHMC) algorithm [llj, which is free of finite step 
size errors, we find that the bare quark mass am^ is reduced by 25%, and the physical 
pion mass by 10%, compared to the accepted values determined previously using the 
R-algorithm. We then extend our simulations to cover a wide range of quark masses, 
mapping out the critical line up to the neighbourhood of the physical point. In agree- 
ment with expectations, the physical point is found to be on the crossover side of the 
boundary. Assuming 0(4) behaviour for the Nf = 2 chiral limit, the fit to our criti- 
cal line can be extrapolated to the m Ut d = axis consistently with the required 0(4) 
scaling behaviour, putting the tri-critical point in that scenario (see Fig. [1]) around 
m*" c /T ~ 2.8. However, non-0(4) behaviour is not excluded by our data. Our results 
should also provide a testing ground and input for analytic attempts to determine the 
critical line from effective theories based on universality arguments (12] (for a review, 
see [13]). 

In a second set of simulations, we repeat the analysis for an imaginary baryon chem- 
ical potential /xb/(zT) = 2.4 and determine the corresponding shift of the critical line, 
following the strategy already used in [9] . Together with additional imaginary \i simula- 
tions for the Nf = 3 case, this allows for a determination of the curvature of the critical 
surface at fiB = 0, which can be readily continued to real values of [in- We find this 
curvature to be negative, as illustrated in Fig. [2] (right). In the (T — //) phase diagram 
this implies that the critical endpoint moves to smaller fi with growing quark mass, 
until it disappears entirely for physical quark masses. This is contrary to customary 
expectations, and in contradiction with the results of [H] , obtained at the same lattice 
spacing and with the same action, but using the R-algorithm and a different numerical 
approach. Clearly, a careful study of systematic errors, due in particular to the very 
coarse lattice spacing, is needed. Still, if the physical point of QCD is indeed in the 
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crossover region at /xg = 0, our finding would imply that the transition will remain an 
analytic crossover also for any finite /ig < 500 MeV, placing a possible QCD critical 
point at much larger values of [ib- Preliminary results to this extent have already been 
given in [15j. 

After summarising the properties of QCD at imaginary \x and introducing the Binder 
cumulant as our observable for the order of the phase transition in Sees. El EJ respec- 
tively, we begin our analysis in Sec. H]with a thorough discussion of step size effects for 
Nf = 3 and a comparison of results from the R- and RHM C-algorithms. The computa- 
tion of the chiral critical line for /ig = and /Jb/ (zT) = 2.4 is presented in Sec. |5l which 
also discusses the resulting new scenario for the (T — hb) phase diagram of physical 
QCD. An assessment of systematic uncertainties is contained in Sec. [HI along with our 
conclusions. 

2 QCD at imaginary \i 

In order to study the phase diagram Fig. [2] at finite baryon density, we employ sim- 
ulations at imaginary chemical potential \x = ifa, where the fermion determinant is 
positive, followed by analytic continuation, as discussed in detail in previous work 
|16[ 15]. To render the paper self-contained, we briefly recall some points needed in the 
sequel. The QCD partition function at finite baryon chemical potential \ib = 3/x is 
even under reflection /i — ■> — /j. Moreover, it is periodic in the imaginary direction, with 
period 2n/N c for N c colours [IT], i.e. Z((j, r /T, }ii/T) = Z(fi r /T, ^i/T + 2tt/3). Because 
of the fermionic boundary conditions, this symmetry implies that a shift in /ij by 27r/3 
is exactly compensated by a ^^-transformation, so that Z(3) transitions take place 
between neighbouring centre sectors for all (/ij/T) c = ^ (n + , n — 0, ±1, ±2, .... It 
has been numerically verified that these transitions are first order for high tempera- 
tures and a smooth crossover for low temperatures [HI [18], as in Fig. [31 Hence, the 
first of these transitions limits the radius of convergence for analytic continuation to 
the first sector for most observables. With a pseudo-critical temperature of T ~ 170 
MeV, our accessible physics range is thus \lb ^,500 MeV. 

Within this first sector, observables can be simulated at imaginary \x = ifa. The 
results may be fitted by truncated Taylor series (O) = J2n c n(Hi/T) 2n , whose con- 
vergence can be tested by inspection. Analytic continuation of successful fits is then 
trivial. 

A remarkable finding of previous work inspecting Taylor series is that, within ^ 500 
MeV, convergence is rapid and screening masses [HI [20] , the pressure [2TJ as well as 
the pseudo-critical temperature [HI [S] are all well described by the leading term ~ ji 2 . 
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Figure 3: Schematic phase diagram for QCD at imaginary chemical potential. The dia- 
gram is periodically repeated for larger values of //j. The heavy lines indicate first order 
transitions, the thin lines crossovers. The nature of the temperature- driven transition 
depends on the parameters or the theory (Nf, quark masses). 

This becomes plausible by noting that at finite T the natural expansion parameter 
is fi/(nT) rather than \ijT [T^JED]- Hence we write the pseudo-critical temperature, 
separating the hadronic from the plasma phase, and the critical quark mass marking 
the boundary between first order and crossover behaviour, as 

l + &i(m M ,m s ) ^y + ... (1) 
l + CiK)(^) 2 + ... (2) 

The choice to treat m u ^ as the independent variable parametrizing the critical line 
in Eq. ([2]) reflects our practical procedure, namely to fix m s and then scan in m u ^ 
because the critical line is a steeper function of the latter. We shall determine the 
coefficients b\, c\ quantitatively and free of step size errors for the Nf = 3 theory, and 
provide the sign of c\ along the whole chiral critical line for Nf = 2 + 1. 

3 Universality and the Binder cumulant 

There are different ways to investigate and exhibit critical behaviour of the theory 
along the critical line, such as finite size scaling (FSS) of susceptibilities, of Lee- Yang 
zeroes or of the Binder cumulant. The Binder cumulant [22] offers various advantages 



To{m Utd ,m s ,fi) 
T (m M , m s ,fj, = 0) 
ra£ |d (m a ,//) 

mC u,d( m s, V = 0) 



for our particular study. It is defined as 

B^m ,/i ) - ,, 5X ^ 2 , (3J 

with the fluctuation 5X = X — (X) of some observable X around its mean value, 
evaluated at the pseudocritical coupling (determined by a peak in ((5X) 2 ) or a zero in 
((5X) 3 )). In practice we study Xi = = 1,2 for the two quark masses m u ^m s . 

In the infinite volume limit the Binder cumulant behaves discontinuously, assuming 
the values 1 in a first order regime, 3 in a crossover regime and some critical value 
reflecting the universality class at a second order critical point. On a finite volume the 
discontinuities are smeared out and flattened, so that B 4 passes continuously through 
the critical value. The location of the critical point in parameter space will then be 
displaced by some finite volume correction. 

In [7] the Binder cumulant was chosen as observable because its critical value for the 
expected Z(2) universality class is distinct from those corresponding to other symme- 
tries like 0(2), 0(4) etc. In this way Z(2) scaling for the Nf = 3 chiral critical point 
was clearly established in [7j. Once the universality class is ascertained, the Binder 
cumulant allows to approximately map out the critical line on a fixed lattice size, by 
scanning the parameter space for the line on which _B 4 is held constant at its critical 
value. In practice this is best achieved by holding one quark mass fixed, and scanning 
in the other. In the small m u ,d regime the critical line turns out to be a steep function 
m s{ m u,d), and thus we choose to scan in the light quark mass while keeping m s fixed. 
In the neighbourhood of a critical point -B 4 can then be expanded in a Taylor series, 

B 4 (am u4 , am s , a/x) = ^ b nm {m s ) (am U)d - am c u d (m s )) n (a/i) m , (4) 

n,m 

with boo(m s ) — > 1.604 for V — > oo. For Nf = 3 there is only one mass variable 
^u,d = n^s — m in the above expression, and the mass dependence of the coefficients 
disappears. 

For large volumes the approach to the thermodynamic limit is governed by univer- 
sality. Near a critical point the correlation length diverges as £ ~ r~ u , where r is 
the distance to the critical point in the plane of temperature and magnetic field-like 
variables, and v « 0.63 for the 3d Ising universality class. In practice, we first find the 
pseudo-critical gauge coupling j3 for a given pair (m Ujd ,m s ), and then compute -B 4 for 
those parameter values. Since (3 is tuned to (3q always, we have r = \am u ^ — am c u d (m s )\. 
_B 4 is a function of the dimensionless ratio or equivalently (L/^) 1 ^. Hence one 
expects the universal scaling behavior 

B, {(L/0 1/u ) = B 4 (L^iam^ - am c u4 (m s )) . (5) 
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4 Nf = 3 without step size errors 



The algorithm most widely used in simulations of finite temperature QCD in the stag- 
gered formulation, both standard and improved, is the R-algorithm [10], which was also 
employed in previous studies of the chiral critical point for Nf = 3 El E] . As pointed 
out in [10], the "correct" usage of the R-algorithm consists of performing simulations 
for various choices of decreasing stepsizes, followed by an extrapolation to zero stepsize. 
However, in practice usually a shortcut avoiding the extrapolation is applied: for some 
reference value of the quark masses, choose a step size for which the step size error is 
smaller than the typical statistical error of the simulation. The molecular dynamics of 
the R-algorithm then suggests to keep the ratio of quark mass am q and step size 5t 
constant, i.e. adjust the stepsize accordingly when the quark mass is reduced. A typical 
choice is 5t = |am g , although in many cases |am g or even am q have been adopted. 

While this procedure has been followed successfully in the intermediate quark mass 
regime, it breaks down for small quark masses, where the linear relation no longer 
appears to hold and the step size needs to be decreased faster than proportionally. 
Furthermore, in a study of the QCD phase transition at finite isospin chemical po- 
tential it has recently been demonstrated that a finite step size leads to a systematic 
underestimate of B4 [23]. Hence too coarse a stepsize can fake a first order transition, 
when the zero step size result really represents a crossover behaviour. 

4.1 The order of the transition: R- vs. RHMC algorithm 

In order to control this important source of systematic error, we have returned to our 
investigation of the Nf = 3 critical point at fiB = [9J, this time with the RHMC- 
algorithm. This algorithm has no stepsize errors and is exact. For a discussion of the 
algorithm and numerical test results, see [24] . Our numerical procedure to compute 
the Binder cumulant is as follows. For each set of fixed quark mass and chemical 
potential, we determine (3q by interpolating from a range of typically 3-4 simulated 
(3- values by Ferrenberg-Swendsen reweighting [25]. For each simulation point 50k-200k 
RHMC trajectories have been accumulated, measuring the gauge action, the Polyakov 
loop and up to four powers of the chiral condensate after each trajectory. Thus, the 
estimate of B4 for one mass value consists of at least 200k, and the estimate of the 
critical mass of at least 800k trajectories. 

Fig. H] (left) shows results of this first study, comparing measurements of B± from 
the RHMC algorithm with those obtained from the R-algorithm at various step sizes. 
The figure confirms the finding from [23J that decreasing the step size leads to an 
increase in the values of the Binder cumulant. It also constitutes a useful test of the 
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Figure 4: Left: Comparison of the Binder cumulant computed with the RHMC algo- 
rithm (leftmost data) and the zero stepsize extrapolation of the R-algorithm. The solid 
lines represent a common fit to all data, the vertical line marks the commonly used 
R-algorithm step size when no extrapolation is performed. Right: Determination of 
m c (fi = 0) = itLq with the RHMC algorithm. The arrow marks the result from the 
R-algorithm JM/. 

RHMC algorithm, whose results indeed correspond to the zero step size limit of the R- 
algorithm. We note that for our smallest quark masses studied, am = 0.005, the RHMC 
algorithm runs over 20 times faster than the R-algorithm at the commonly applied step 
size. Since the latter also requires runs at several step sizes for the extrapolation, the 
RHMC is thus considerably more economical in producing results free of step size 
errors. 

4.2 The Nf = 3 critical quark mass at fi = 

In order to eliminate step size errors, we now proceed to repeat the calculation of the 
critical quark mass mjj in the three-flavour theory by means of the Binder cumulant, 
this time with the RHMC algorithm. The result obtained on an 8 3 lattice is shown in 
Fig. H] (right). Qualitatively the behaviour is the same as previously, with B4 growing 
from first order behaviour through its critical Ising value to crossover with increasing 
quark mass, which can be fitted to leading order in the quark mass. However, the 
critical Ising value is now obtained at a bare mass of amjj = 0.0260(5), which is about 
25% smaller than the value am^ ph 0.033 quoted by all previous work using the R- 
algorithm 

One may ask whether this change affects bare quantities only, while the R and 
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Figure 5: Finite size scaling of the Binder cumulant for Nf = 3. The lines represent a 
common fit to the data, according to Eq. (Ti|). 

RHMC algorithms probe the same physics. To study this issue, we measured the 
zero-temperature hadron spectrum at the parameters (f3 C) am c ) RHMC using RHMC, 
and compared with the same exercise performed in [7] at the parameters (/3 C , am c )# 
using the R-algorithm. For the pion, which is the most accurately determined, the ratio 
mjTo changes from 1.853(1) (R [7J) to 1.680(3) (RHMC). This reduction of 10% in the 
pion mass corresponds to a change of 20% in the renormalized quark mass, very near 
the observed 25% change in the bare quark mass. Therefore, replacing the R by the 
RHMC algorithm corrects a large error in the physical values of the critical parameters. 
The correction should be even larger for smaller m u ^ masses. We conclude that for 
the study of the QCD phase transition in the region of physical quark masses, step 
size errors in the Monte Carlo algorithm can lead to a qualitatively different picture at 
fixed parameter values, and the use of an exact simulation algorithm is mandatory. 

Our results so far have been obtained on a single spatial volume 8 3 . The next 
task is to study the FSS behaviour and uncover possible finite volume effects. This 
is particularly important since large finite size corrections were reported in a recent 
investigation of the chiral critical point at finite density in the Taylor expansion of 
susceptibilities (26]. Fig. [5]shows data obtained for three lattice sizes with L — 8, 12, 16. 
According to Eq. §5§ and the corresponding discussion, in the scaling region near a 
critical point B± should be described by 



B 4 (m,L) = b + bL 1/v ( 




(6) 



9 



5.12 



-0.01 







0.01 0.02 0.03 0.04 0.05 0.06 0.07 



(aw) 2 



Figure 6: The pseudo-critical coupling /3 Q for Nf = 3 on an 8 3 x 4 lattice as a function 
of quark mass and imaginary chemical potential. For each value of fj, i; f3 is determined 
from 4 data points corresponding to masses am = 0.02 — 0.034. The two lines represent 
the fits from TableUl 

with bo = 1.604 and v = 0.63 for 3d Ising universality. We have checked for finite 
volume effects by fixing 6 to the Ising value and fitting for b, v and m$. With a \ 2 of 
0.74 per d.o.f., we obtain am^ = 0.0263(3) consistent with our result from L = 8 only, 
and v = 0.67(13), which is consistent with the Ising exponent. We conclude that for 
Nf = 3 the Binder cumulant is close to thermodynamic scaling for lattice sizes L > 8, 
and hence finite volume effects are under control in this calculation. 

4.3 The pseudo-critical temperature for Nt = 3 at finite fi 

On the lattice, T is determined from the pseudo-critical gauge coupling, which we 
define as the location of the peak of the plaquette susceptibility. On any finite volume 
it can be expanded as a double series in mass and chemical potential around the three 
flavour critical point m^, 



Fig. [6] shows the measured values of /?o(am, a/Xj) for four different imaginary chemical 
potentials spanning the whole /ij range up to the first Z(3)-transition. For each value of 
a fii, four quark masses in the range 0.02 < am < 0.34 have been simulated. As in our 
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Table 1: Fzts of the Taylor expansion (3o(am,an), Eq. ([?]), cmr data in Fig. 



previous study, we obtain good fits retaining the leading terms only, as shown in Tabled! 
In particular, the term quadratic in the chemical potential is now sufficient to describe 
the data all the way out to fi/T = 1, the quartic coefficient being consistent with zero. 
Using the two-loop beta function, this translates to a pseudo-critical temperature of 

T °<"^ - 1 + 2 . 111(17 )f!!i^»Vo.667(6)f^) 2 + 0.23( 9 )(^) 4 + ... (8) 



'/u("'o-0) ~ V kT 7 V ; ' ;/ n/ V-'/j. 

Again, we note a shift of up to 10% in the coefficients compared to the R algorithm 
results [9]. Eq. ([8]) has to be considered with some caution, since it is well known 
that the two-loop beta function is rather inaccurate at our coarse lattice spacing. The 
effect of the non-perturbative beta function is to increase the absolute values of the 
coefficients A, B, perhaps by up to a factor of 2 [2 



4.4 Quark mass dependence of the critical point for Nf = 3 

In order to detemine the order of the transition, we now repeat the previous procedure 
to find how the critical bare quark mass m c (n) changes with imaginary chemical po- 
tential. As in the case of the pseudo-critical temperature, we express this by a Taylor 
series Eq. ([2]): 

^0)— (&)' + -. < 9 > 
in order to be able to continue to real fi. Inversion of this function will then give the 
location of the critical point as function of the quark mass, /i c (m). In practice, at a 
finite lattice spacing we are dealing with the expansion in lattice units, 

am c (ajj,) = am c (ajj, = 0) + c'^a/i) 2 + c' 2 (a/i) 4 + . . . (10) 

A crucial point is that, for fixed temporal lattice extent N t , the lattice spacing entering 
the dimensionless am c (/i) and am c (0) is different, since T (m c (//), //) = l/(N t a(fj,)) 
depends on fi. The relation of the leading coefficient to its continuum counterpart is 
thus given by 

1 dm c 7r 2 c[ 1 G?To(?n c (/i), /u) 

Cl = ^ rf(/i/vrT) 2 = ivf ^ + T Q {m c , 0) rf(/i/vrT) 2 ' ^ 
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or in terms of A = 2.111(17), B = -0.667(6) from Eq. © and d x from Eq. (HDD, 



^l^^l 1 -^' ' (12) 



'o 



The coefficients are extracted from our data for B^ obtained at imaginary ll — ifii, 
by fitting to a double expansion about the known critical point at m c (/i = 0), 

B 4 (am, a fx) = b ni (am — am C Q) n (aLi) 21 . (13) 

n,l 

The leading coefficients c[ are then obtained as 



c , _ c?am c dB A . ^ 4 , _ 

1 d(aLi) 2 d(aLi) 2 \dam J b w 



d 2 am c d 2 B d ( dB A \' x dB 4 ( dB^y 2 d 2 B 4 





(dB 4 
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(dB 4 




{dam 


601611 


ft 2 

°10 



+ 



d[(aLi) 2 ] 2 d(aLi) 2 \dam J d(aLt) 2 \dam J d(aLi) 2 dam 

602 . 6 i6ii 

= — b — + ^r- ( 15 ) 
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For the actual analysis it is thus convenient to reparametrise the second order expansion 
of -B4 as 

5 4 (am, a/i) = 1.604 + b w [am — am$ — c'^a/z) 2 ] + 6 2 o( am — am^) 2 

—610 [(c' 2 — d 1 C)(aLi) 4L + C(am — am^^aLi) 2 ^ , (16) 

with C = —611/610, and fit the data via the parameters mjj, b w , b 2 o, c[, c' 2 , C. 

Our data for five different values of imaginary chemical potential are shown in Fig. [7] 
(left). Remarkably, there seems to be negligible influence of the chemical potential. 
The results of various simultaneous fits of all four curves are displayed in Table El 
All fits are good, and none of the next-to-leading terms is significantly constrained. 
This is corroborated by discarding all next-to-leading terms, which leads to a perfectly 
acceptable fit with c\ consistent with zero, as in the last line of Table El Fig. [7] (right) 
displays the error band coming from a linear fit (Table El line 3). Clearly, the slope c[ 
is very nearly zero. 

The final result is then obtained by employing Eq. (fl~2|) to convert to continuum 
units. The second factor in Eq. (fl~2l) is 1.077(2), close to 1. In the first factor, the term 
B, which describes the variation of Tq(ll) with real 11 and is thus negative, reinforces 
the negative trend of c[, to yield 
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Figure 7: Left: B^[am^a[x) for different imaginary chemical potentials. The lines 
correspond to the simultaneous fit of all data according to Eq. [To]) and Table (H line 
3. Right: One sigma error band from a linear fit (Tabled line 3) of the critical quark 
mass as a function of imaginary \x. 



m 


b w 


&20 




4 


C 


X 2 /dof 


0.0262(7) 


13.3(1.4) 


-91.6(143.5) 


-0.079(47) 


-1.6(1.0) 


-2.1(3.5) 


0.90 


0.0263(6) 


13.9(0.6) 




-0.075(42) 


-1.35(0.73) 




0.82 


0.0270(5) 


13.6(0.6) 




-0.0024(160) 






0.93 


0.0271(3) 


13.6(0.6) 










0.88 



Table 2: Fitting B^am, afi) to a Taylor expansion to different orders in the indepen- 
dent variables, according to Eq. [To]) . The numbers of d.o.f. are 14,16,17,18, respec- 
tively. 



Hence, we arrive at the surprising result that the first order region in the phase diagram 
Fig. [2] shrinks when a real chemical potential is switched on. This is contrary to the 
expected qualitative behaviour. 

The reader will notice the large error on the coefficient in Eq. ( fTTj) . It is a conservative 
estimate and stems entirely from the larger error on c[. If one were to include a /i 4 - 
term, the previous conclusion would only be strengthened: the leading term gets more 
negative and a negative quartic term comes on top of it, cf. Table |2j 

However surprising, our findings agree with preliminary results for the same lattice 
theory at finite isospin chemical potential, which indicate that there too, the transition 
becomes weaker as the chemical potential is turned on {28] . Finally, let us note that the 
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same qualitative behaviour applies to the first order region in the heavy quark limit, 
which has recently been shown to also shrink with real chemical potential [29] . 

5 The chiral critical surface for Nf = 2 + 1 

Having removed finite step size errors from the Nf = 3 calculations, we proceed to map 
out the chiral critical line for non-degenerate quark masses. All our simulations have 
been performed with the RHMC algorithm. Since 8 3 x 4 lattices proved to be large 
enough for our observable in the case of Nf = 3, we use that lattice size to trace out 
the critical line, performing another check of finite volume effects at am u ^ = 0.015 on 
a 12 3 x 4 lattice. 

With two different quark masses in the theory, a technical question concerning the 
Binder cumulant arises. Obviously, B± can be constructed from the chiral condensate 
of either mass flavour. Universality guarantees that, in the infinite volume limit, either 
choice tends to the same universal value. However, in a finite volume there are cor- 
rections, and they are different for different operators. The corrections are minimised 
for that superposition of operators, which corresponds most closely to the mapping of 
the QCD parameters onto the scaling fields of the effective 3d Ising model. It is well 
known that, even for the case of three degenerate flavours, this is a superposition of 
the chiral condensate and the plaquette [7j, as well as higher dimension fermionic and 
gauge condensates. 

Here, we do not attempt to construct an optimised observable by mixing in gauge 
condensates, but simply compare the behaviour of _B 4 constructed from condensates of 
different mass flavours, Xi = ipiipi, as shown in Fig. El The observables B±(X S ), Bi(X u ^) 
constructed from the condensates of the heavy and light flavours, respectively, are ob- 
served to intersect the critical value at significantly different values of am Uj d- Neverthe- 
less, comparison of results obtained at L — 8, 12 shows that this difference is rapidly 
disappearing on larger volumes. Moreover, the common intersection point to which 
the results converge appears to be close to that obtained from _B 4 (X Ujrf ) on L = 8, 
indicating that the latter has far smaller finite volume corrections. This is not too 
surprising, as one would expect the scaling field corresponding to chiral symmetry to 
be dominated by the lightest quark flavour. Hence, in the following we will always 
work with 5 4 (X M ^) constructed from light quark condensates. 
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Figure 8: Finite size scaling of the fits to B4 at fixed am s = 0.075, constructed from 
the light (solid line) and the heavy (dotted line) flavour condensates for L = 8, 12. The 
light flavour observable has much smaller finite volume corrections. 

5.1 The critical line for Hb = 

By fixing m s and scanning in m u ^ (at least 4 values), the critical light quark mass 
for that choice of m s is determined by interpolation, analogously to the three-flavour 
case. This is repeated for other values of m s , resulting in the sequence of critical points 
m ud( m s) displayed in Fig. left. As in the three-flavour case, every critical point 
appearing in this figure consists of at least 800k RHMC trajectories. 

There are several features of Fig. worth discussing. An interesting observation 
concerns the behaviour of the function m c s (m u ^) in the neighbourhood of the three 
flavour critical point. If is constructed from gauge condensates and neglecting the 
change in the pseudo-critical coupling Po( m u,d, m s) with the quark masses, a Taylor 
expansion around the symmetric critical mass t?1q yields for the line of constant (critical) 
B4 the leading order result [7j 

m s = ml-2(m U)d -m c fj ), (18) 

i.e. the critical line should pass through the symmetrical point with slope -2. In con- 
trast, our data extracted from B^Xu d) exhibit a different slope, see Fig. [9] (left). This 
underlines again the importance of choosing an appropriate observable for finite vol- 
ume computations. A Taylor expansion of B4 as in Eq. f[T8l is only possible on finite 
volume, but expanding B^(X U ^) would yield additional non-perturbative contributions 
to Eq. (|T8|) . We thus conclude that Eq. ({TBI does not describe the critical line, not 
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am ud am ud 



Figure 9: Left: The chiral critical line in the bare quark mass plane at fig = 0. The 
heavy line indicates the Nf = 3 diagonal. Also shown is the physical point according 
to llJ$ , and a fit to extrapolate the line to a possible tricritical point on the m s -axis. 
The arrows mark the points where T = simulations were performed to set the scale, 
Sec. I5.H Right: Comparison of the critical line at \Lb — and /^/(zT) = 2.4. 

even in the immediate neighbourhood of Nf = 3. 

Another interesting question is how the critical line continues to even smaller light 
quark masses. If the chiral limit of the Nf = 2 theory exhibits 0(4) universality, 
then the critical line hits the axis m Uj d = in a tricritical point at some finite strange 
quark mass value m*™ c [2]. Whether this scenario is realized or not is an issue not yet 
settled (cf. the discussion and references in [H]). Among the most recent publications 
using staggered fermions, one favors a first order scenario for the chiral limit [30J while 
the other supports the 0(4) scenario [3T] . With our current data, we are unable to 
decide this question, but we can check for consistency with the 0(4) scenario, which 
implies mean-field exponents near the tricritical point (m u ^ = 0,m s = m^" c ). Indeed 
our data support a fit to the ~ m u d approach to the chiral limit, as shown in Fig. [9j 
left, predicting the tricritical point to be at am* riC ~ 0.7 or m*" c /T ~ 2.8. Note 
however that (i) our N t = 4 lattice is very coarse (a ~ 0.3 fm), and (ii) our spatial 
volume becomes rather small as m u ^ is reduced: for the uppermost point in Fig. [9], 
left, corresponding to the physical strange quark mass, m^L ~ 1.7 only. Thus, our 
systematic error might be rather large. Nevertheless, we have strong indications that 
mf lc is significantly larger than the physical strange quark mass. 
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(am u , d , am s ) (3 


am n arriK am p 




(0.0265,0.0265) 5.1374 
(0.005,0.25) 5.1857 


0.420(1) 0.420(1) 1.383(7) 
0.2109(1) 0.8915(1) 1.398(16) 


0.304(2) 0.304(2) 
0.151(2) 0.638(8) 



Table 3: Parameters and meson masses for the low-temperature (12 3 x24) simulations 
performed to set the scale. The mass ratios of the second line imply that this point 
corresponds to m s = {m s ) phys ,m u>d < (m U)d ) phys . 



5.2 The chiral critical line and the physical point 

The most important question regarding the critical line is, of course, its location relative 
to the physical point of QCD. So far, all known lattice data are "consistent" with 
the physical point being in the crossover region. This is also the result found by a 
simulation at physical quark masses in [H]. However, these results were obtained by 
the R-algorithm, and we have seen in the three flavour case that significant shifts in 
the critical quark masses can arise due to step size errors. 

In order to estimate the location of our critical line in physical units, we have there- 
fore performed zero temperature simulations with bare quark masses corresponding to 
two points on the critical line. One corresponds to the three flavour theory and the 
other to the point with roughly physical strange quark mass, as indicated by the arrows 
in Fig. [9] (left). The parameters of the simulations are given in Table [3j together with 
the measured meson masses. In both cases, the lattice size was 12 3 x 24, and about 
400 configurations were analyzed. 

Setting a physical scale along the critical line is a tricky problem. Neither of our 
simulation points matches the physical (m u> d,m s ) point, so that strictly speaking one 
cannot match to any real world observable. Doing so anyway, different observables 
inevitably give different values for the lattice spacing. A measurement of the qq force 
via elongated Wilson loops gives r /a = 1.85(2) and 1.87(2) respectively, where r = 0.5 
fm is the Sommer scale. This amounts to a(ro) ~ 0.27 fm in both cases. On the other 
hand, matching the p-mass to its physical value gives a lattice spacing which is by 20% 
larger. Note, however, that a difference of similar magnitude has been observed in [TJj 
on the physical point. This suggests that the greater part of this difference is due to 
cut-off effects rather than to the deviation from physical parameters. 

It thus appears safer to avoid setting an absolute scale altogether, and instead com- 
pare the meson mass ratios from Table [3] with their values at the physical point, 
(m 7T /m p )phys = 0.18 and (mK/m p ) p h ys = 0.645. We thus conclude that our Nf = 2 + 1 
point on the critical line indeed corresponds to the physical Kaon mass and to pions 
lighter than physical. In other words, the physical point is on the crossover side of the 
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Figure 10: Comparison of our Nf = 2+1 meson masses aM (leftmost points) with 
those obtained in Ref. \Tl$ with the R-algorithm. Good consistency is apparent. 



critical line. 

It is interesting to compare our Nf = 2 + 1 meson masses, in lattice units, with 
those of [H] , which used the same strange quark mass, and the R-algorithm at almost 
the same inverse coupling (3 = 5.19. The comparison is shown in Fig. [101 where the 
straight lines are fits to the data of [2] only. Good consistency is apparent, showing 
that the R-algorithm step size error for the meson masses is small, for the parameters 
considered in |14J §. The figure also shows mK/m p to be practically independent of 
m Uj d, thus affirming our conclusion above. 

Finally, the fact that the lattice spacing varies little between our two simulation 
points, implies that T c itself does not change much as one moves along the critical 
(m n rf, m s ) line. This is in agreement with model calculations [T3] . 

5.3 The critical line for /ig/(zT) = 2.4 and the critical surface 

We have also run a second set of simulations with an imaginary chemical potential 
HB/(iT) = 2.4, in order to determine how the critical line shifts with baryon density. 
These data are shown in Fig. (right), in lattice units. We observe that the shift in the 
critical line is very small, despite the sizeable value of the chemical potential. Within 
two sigma the line is consistent with its //# — counterpart. Moreover, to the extent 

2 Note also, that there are preliminary results by Z. Fodor and S. Katz 
(https://www.bnl.gov/sewm/) using finer lattices and the exact RHMC algorithm, which con- 
firm that the physical quark masses give a finite-temperature crossover. 
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coo = Po(0,m c u>d ) 


cio, (/i 2 ) 


C20, (/^ 4 ) 


coi, (m) 


X 2 /dof 


5.1838(3) 


0.572(9) 




1.75(13) 


1.5 



Table 4: Fzt 0/ t/ie Taylor expansion /3o(am UtC [, an) , Eq. ([?]), to our data for fixed 
am s = 0.25 and am c ad ~ 0.0064. 

that there is a displacement, it shows a trend to lie to the right of the zero density line. 
This qualitative observation is in accord with our earlier finding in the three flavour 
case (Eq. [10]) that c[ ~ or slightly negative: in lattice units, the first order region 
tends to expand slightly as an imaginary chemical potential is turned on (see Fig. [7] 
(right)). 

Similarly to the Nf = 3 case, one expects the data along the whole line to be 
described by the leading term in the Taylor expansion, 

am c u d (am s , a/i) = am c u d {am s , a/i = 0) + c'^arris) (a/i) 2 + . . . (19) 

where now c[ depends on am s . Conversion to continuum units proceeds as for Nf = 3, 
by determining the equivalent of Eqs. (l8"][T2l for Nf = 2 + 1. We do this for fixed 
physical strange quark mass am s = 0.25 and scanning in m u ^ which now plays the 
role of the variable quark mass, and find om c ud ~ 0.0064. For the variation of the 
pseudo-critical coupling, Eq. ([7]), we obtain the coefficients given in Table HI leading 
to A = 1.90(13) and B = —0.49(1) in this case. These are similar in magnitude to 
the three-flavour values. Note that c[ ~ means that m° u d /T c remains constant as the 
chemical potential is turned on. The decrease of the pseudo-critical temperature with 
real fi, given by B, is then the dominant effect. It dictates that m c ud decreases when a 
real chemical potential is turned on. In other words, the first order region shrinks. 

While our data for Nf = 2 + 1 and small quark masses have larger errors which do 
not yet allow to constrain Ci{m u ^) quantitatively, Fig. [9] and Eq. ({T21 leave little doubt 
that this coefficient is going to be negative along the whole upper part of the line. We 
thus arrive at the conclusion that, for the lattice spacing considered here, the curvature 
of the critical surface at fiB = is negative, and the first order region is shrinking when 
a real chemical potential switched on. 

5.4 An alternative scenario for the QCD phase diagram 

Let us take our results at face value for a moment and consider the implications if such 
a qualitative result also holds in the continuum limit. This leads to a scenario for the 
(T, /i) phase diagram which is at odds with common expectations. We find that the 
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Figure 11: For dm c (n)/dn 2 < 0, there is no critical point at all, the dotted line on the 
right is merely a crossover. Any additional critical structure would not be continuously 
connected to that at fi = 0. 

first order region in a plane of constant /is is actually shrinking with growing real \ib- 
If the physical point is in the crossover region at \ib — 0, then switching on a chemical 
potential will not lead to an intersection with the critical surface as long as the latter 
is well described by its curvature at \i& = 0, i.e. for fi/T < 1, or equivalently \ib < 500 
MeV. In the absence of any additional (and so far unknown) critical structure, there 
would then be no critical point or first order phase transition at all. The (T, //) phase 
diagram of physical QCD would then only have the possible transition line separating 
the superconducting phase from nuclear matter, as in Fig. [11] (right). 

Note that this scenario is perfectly consistent with all universality arguments and 
the known results for \x = 0. This can be illustrated in the three flavour theory by 
considering the change of the (T, /i)-diagram with quark mass, as depicted in Fig. O 
All boundary conditions are met, in particular there is a first order phase transition at 
H = for quark masses smaller than m§. However, according to the negative sign for 
Ci in Eqs. f|T7ll9l) . the critical endpoint is now moving to the left with growing quark 
masses, until it disappears entirely. 

It is natural to ask how reliable this unexpected scenario is regarding systematic 
errors. We have seen in Sec. 15.21 that N t = 4 lattices are very coarse, and we have dis- 
cussed the enormous sensitivity of the critical values of the mass parameters (m U) d, m s ) 
to cut-off effects [32j. It is therefore expected that the values for m c (/Xi) shift once this 
analysis is repeated on finer lattices and/or with improved actions. The crucial ques- 
tion is whether the continuum extrapolated slope c' x in Fig. fright) will be positive or 
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Figure 12: The Nf = 3 phase diagram as a function of quark mass. The critical 
endpoint now is moving to smaller \i c with growing quark mass, due to the negative 
sign of c\ in Eqs. ( f!7|fgj) . 

negative and how it will balance out against the contribution from the pseudo-critical 
temperature. Or, equivalently, how the ordering of the zero and finite \x critical lines 
in Fig. [9] turns out in the continuum limit once physical units are used. In particular, 
our continuum conversion is sensitive to the non-perturbative beta function. A look 
at Eq. (TT2"|) gives a quick estimate of what is needed for a positive C\. Either c[ would 
have to be positive of the order ~ am^Nljv: 2 , or A has to grow by a factor larger than 
10 on the way to the continuum. Thus, while we presently cannot rule out that the 
picture reverts back to the standard scenario in the continuum limit, the opposite is 
obviously also possible as suggested by our data. 

Finally, all of our arguments are based on the simplest scenario, in which the finite 
H critical point of physical QCD is continuously connected to a critical point for some 
other mass values at /x = 0. We cannot exclude a more complicated possibility, where 
the phase boundary of the color superconducting phase (see Fig. [12]) would distort as a 
function of the quark masses, and give birth to a critical point distinct from the one we 
study, and which would survive for physical quark masses. This would correspond to 
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a scenario with an additional critical surface in Fig{TT] (left) above the one we studied 
here. We thus conclude that even the qualitative features of the QCD phase diagram 
cannot be regarded as settled yet. 

5.5 The finely tuned critical end point 

While based on our present data we are unable to make a reliable quantitative predic- 
tion for the location of the critical point for physical QCD, we have obtained important 
qualitative information regarding its quark mass dependence. Irrespective of the con- 
tinuum extrapolated sign, all our evidence is that the absolute value of c\ in Eq. (J2J) 
will be ~ 0(1) as naturally expected, while the effect of subleading terms is small up 
to hb ~ 500 MeV. This means that the critical quark masses are very weakly varying 
functions of the chemical potential, which is in line with the corresponding behaviour 
of the pseudo-critical temperature or, indeed, the equation of state. Consequently the 
inverse function // c (m U) d) is very strongly varying with quark mass, i.e. the finite \i 
critical surface emerges very steeply from the \x = critical line in Fig. [2j Hence, even 
if in the continuum limit the conventional scenario with positive curvature is realized, 
the precise location of the critical end point will be exceedingly quark mass sensitive. A 
simple estimate using c± = 1 shows that, in order to have [i° B < 400 MeV, the physical 
point has to lie within < 5% of the chiral critical line, i.e. the physical quark masses 
would be fine tuned. While there is nothing forbidding such a situation, it appears 
rather unnatural. Moreover, if realized in nature, it would make a quantitative de- 
termination of n B through simulations exceedingly difficult. (For example, one might 
even need to treat the u- and <i-quarks as non-degenerate). 

6 Discussion and conclusions 

As we have demonstrated in the preceding sections, step size errors are eliminated 
from our calculations. Finite volume effects are under control in the Nf = 3 theory 
and for am Uj d > 0.015 in Nf = 2 + 1, where we have performed our finite volume 
check. Note that this corresponds to pion masses larger than physical. For the point 
on the critical line with the lightest quark mass considered, we only have m^L ~ 1.7, 
and finite volume effects are to be expected. Ideally the part of the critical line in the 
neighbourhood of the physical point should be checked on a larger lattice. 

However, by far the largest source of uncertainty is due to the coarse lattice spacing 
a ~ 0.3 fm, as evidenced by several aspects of these computations. Strong cut-off effects 
reveal themselves when attempting to set a physical scale for the problem. Moreover, 
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a change in the discretization of the Dirac operator on a lattice this coarse can change 
the pion mass corresponding to the second-order transition by a factor ~ 4 [32] at 
Hb = 0. Finally, it has recently been pointed out that staggered simulations at finite 
Hb suffer from additional discretization errors compared to hb = [33] when Nf ^ 4, 
due to the eigenvalue structure when taking the fourth root of the determinant. For 
simulations at imaginary fig, the eigenvalues are pure imaginary, and this additional 
error is of 0(a 2 ), with possibly a large coefficient. A safe strategy thus is to first take 
the continuum limit of imaginary simulations, and then continue to real fis- For 
reweighting approaches at real fis one even expects 0(a) errors. 

In interpreting our findings and comparing with other work, it is important to take 
systematic uncertainties into account. Given the cut-off effects, the sensitivity of the 
critical point to step size errors and, most notably, to the quark mass, it is clear 
that the discrepancy between our findings and those of [H] is nothing remarkably 
unusual, but merely reflects the large and different systematic uncertainties afflicting 
these calculations. In particular, in [T4] the quark masses am q were held fixed in 
lattice units while afi was increased. Equivalently, m ? /T was kept fixed. However, T 
decreases under the influence of a chemical potential, in a manner similar to Eq. [HJ so 
that the quark masses at the critical endpoint in |H] are about 5-10% smaller than 
physical. This small deviation from a line of constant physics has a large impact 
on the location of the critical point, because of the high sensitivity of the latter on 
quark masses^]. The effect is to artificially move the critical point to smaller chemical 
potentials. The shifted masses may even reside in the first order region, causing a 
critical point to be found even if in fact there is none, consistently with the scenario 
discussed here. 

Our study of the chiral critical surface also suggests that one cannot draw conclusions 
for physical QCD from simulations of the critical point in Nf = 2 QCD [261 EI]- This 
becomes clear when considering Fig. dj which describes the fi = expectations. If 
one moves from the physical point upwards by increasing m s to infinity, the distance 
to the critical line increases considerably. Given the high sensitivity of the critical 
chemical potential /z c to this distance (i.e. the small curvature of the chiral critical 
surface in Fig. [2]) which we observe, one should expect large differences for [f between 
the N f = 2 + 1 and Nf = 2 theories 0. 

Resolving these various systematic issues and deciding which scenario for the (T — fi) 

3 In our approach, the conversion from lattice units Eq. (flQ|) to physical units Eq. Q changes the 
coefficient c[, which is nearly zero, to c±, which is negative. 

4 Moreover, for Nf = 2 one expects a tricritical point n tn (m u( i = 0), with mean-field critical 
exponents which govern the analytic form of the critical line fi(m uc [,m s — oo), in contrast to the 
Nf = 2 + 1 case. 
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phase diagram is realized in nature thus urgently requires further investigations of the 
Nf = 2 + 1 theory with exact algorithms on finer lattices. Among the various finite 
fi approaches, our imaginary \x simulations require comparatively moderate computer 
resources to achieve this goal. 
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